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A  three-dimensional  mathematical  model  coupling  the  electrochemical  kinetics  with  fluid  dynamics  is 
developed  to  simulate  the  heat  and  mass  transfer  in  the  one-cell  stack  of  planar  solid  oxide  fuel  cells 
(SOFCs).  Based  on  flow  uniformity  analysis,  the  distributions  of  temperature,  current  density,  overpo¬ 
tential  loss  and  other  performance  parameters  in  various  operating  parameters  are  obtained  using  a 
commercial  CFD  code  (Fluent)  coupled  with  the  external  subroutines  programmed  by  VC++.  Numerical 
flow  data  are  observed  in  good  agreement  with  experimental  results  reported  in  the  literature.  Results 
show  that  the  one-cell  stack  in  counter  flow  case  has  the  advantages  in  better  uniform  current  density 
and  temperature  distributions  of  PEN  (Positive/Electrolyte/Negative)  structure  in  the  width  direction, 
higher  power  output,  fuel  utilization  factor  and  fuel  efficiency  than  that  in  co-flow  case.  For  counter  flow 
case,  better  thermoelectric  characteristics  are  observed  in  the  temperature  gradient,  power  output,  fuel 
utilization  factor  and  fuel  efficiency  with  the  decrease  in  the  fuel  inlet  flow  rate  or  the  anode  porosity. 
Increasing  the  air  inlet  flow  rate  and  decreasing  the  fuel  inlet  temperature  will  reduce  the  temperature 
gradient;  power  output,  fuel  utilization  factor  and  fuel  efficiency  are  enhanced  with  the  increase  of  the 
air  inlet  temperature  and  the  decrease  of  the  anode  pore  size  and  thickness. 

©  2009  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Solid  oxide  fuel  cell  (SOFC)  has  been  considered  as  a  promis¬ 
ing  alternative  energy  source  device  for  residential  and  distributed 
power  plants  because  of  its  high  energy  conversion  efficiency  and 
power  density,  low  environmental  hazards  and  potentially  low  pro¬ 
duction  cost  [1,2].  However,  the  further  development  of  planar 
SOFCs  faces  the  challenges  related  to  maximize  the  power  density 
and  minimize  the  non-uniform  temperature  distribution,  which 
contributes  to  the  thermal  stress  in  different  SOFC  components 
[3-5].  To  optimize  SOFC  stack  design  and  understand  its  complex 
operating  mechanism,  many  challenges,  such  as  heat  and  mass 
transfer  together  with  electrochemical  reactions,  optimization  of 
geometry  and  development  of  new  materials,  still  remain  to  be 
solved  step  by  step.  The  purpose  of  the  optimal  geometric  design 
is  also  to  get  uniform  distributions  of  air  and  fuel  flows  in  order  to 
improve  the  cell  stack  performance  [6]. 


*  Corresponding  author  at:  State  Key  Laboratory  of  Materials  Processing  and  Die  & 
Mould  Technology,  Huazhong  University  of  Science  and  Technology,  No.  1037,  Luoyu 
Road,  Wuhan  430074,  PR  China.  Tel.:  +86  27  8754  3493;  fax:  +86  27  8754  3493. 
E-mail  address:  xiatianhust@hotmail.com  (W.  Xia). 

0378-7753 /$  -  see  front  matter  ©  2009  Elsevier  B.V.  All  rights  reserved. 
doi:10.1016/j.jpowsour.2009.06.009 


In  the  past,  modelling  the  planar  SOFC  stack  to  study  influences 
of  various  geometric  and  operating  parameters  on  the  stack  per¬ 
formance  has  been  carried  out  under  the  assumption  of  uniform 
flow  velocity  distributions  in  gas  channels  [7-11].  However,  few 
investigations  concerning  the  uniformity  of  air  and  fuel  flows  for 
U-type  and  Z-type  SOFC  stack  have  been  reported  [12,13].  Further¬ 
more,  quit  a  few  researchers  have  also  made  an  attempt  to  study 
the  complex  thermo-fluid  electrochemical  transport  phenomena 
in  planar  SOFC  stack  with  non-uniform  gas  flow  rate  in  the  stack 
direction  and  various  interconnects  [6,14]. 

Many  researchers  have  discussed  the  influences  of  various 
operating  parameters  on  SOFC  performance,  such  as  air  and  fuel 
flow  rates,  anode  thickness,  steam  to  carbon  ratio,  specific  area, 
pre-reforming,  and  others  [15-26].  Therefore,  it  is  necessary  to 
investigate  effects  of  various  geometric  and  operating  parameters 
on  the  distributions  of  air  and  fuel  flows,  such  as  flow  configura¬ 
tion,  delivery  rates  of  air  and  fuel,  inlet  temperatures  of  fuel  and 
air,  and  anode  thickness,  pore  size  and  porosity.  Then,  effects  of  the 
flow  uniformity  of  the  obtained  air  and  fuel  flow  rates  on  the  stack 
performance  are  simulated  to  predict  the  temperature  and  current 
density  distributions,  average  temperature,  average  current  den¬ 
sity,  air  ratio,  power  output,  fuel  utilization  factor  and  fuel  efficiency 
during  the  operating  process  of  SOFC  stack. 
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Nomenclature 

AH  air  channel 

ak  polynomial  coefficients 

B  inertial  coefficient 

Cp  gas  mixture’s  specific  heat  capacity  (J  kg-1  K-1 ) 

CP  j  specific  heat  capacity  of  species  i  (J  kg-1  K-1 ) 

dp  average  pore  size  of  the  electrodes  (pun) 

Dh  hydraulic  diameter  (m) 

Dj  eff  effective  diffusion  coefficient  of  species  i  (m2  s-1 ) 
Di  m  molecular  diffusion  coefficient  of  species  i  (m2  s-1 ) 
DKi  Knudsen  diffusion  coefficient  of  species  i  (m2  s-1 ) 

E  electromotive  force  (V) 

F  Faraday  constant  (C  mol-1 ) 

FH  fuel  channel 

i  local  current  density  (A  m-2 ) 

i0  exchange  current  density  (A  nrr 2 ) 

i0) a  anode  exchange  current  density  (A  m-2 ) 

i0)C  cathode  exchange  current  density  (A  m-2 ) 

Ji  mass  flux  of  species  i  (kg  m-3  s-1 ) 

k  reaction  rate  constant  (mol  m-3  Pa-2  s-1 ) 

I<  permeability  coefficient  ( m2 ) 

L  cell  stack  length  (mm) 

Mi  molecular  weight  of  species  i  (kg  mol-1 ) 

Mm  average  molecular  weight  (kg  mol-1 ) 

na{r  air  feed  header  inlet  flow  rate  (mol  h-1 ) 
ne  electrons  transferred  per  reaction 

nfuei  fuel  feed  header  flow  rate  (mol  h-1 ) 
p  pressure  (Pa) 

Pi  partial  pressure  of  species  i  (Pa) 

Psofc  power  density  of  the  cell  ( W  m-2 ) 

Q  lower  heating  value 

R  universal  gas  constant  (J  mol-1  K-1 ) 

Rj  ohmic  resistivity  ( £2  m) 

Re  Reynolds  number 

S  source  term;  effective  temperature  (K) 

Si  species  source 

Sm  mass  source 

•^T-thermo  energy  source  caused  by  thermodynamic  reaction 
^T-ohmic  energy  source  caused  by  ohmic  resistivity 
Sv  momentum  source 

T  temperature  (I<) 

TaVe  average  temperature  of  PEN  structure  (I<) 
u  velocity  of  rib-channel  (m  s-1 ) 

Ufuel  fuel  utilization  factor 

V  velocity  vector  (ms-1) 

vk  velocity  components  in  k  direction  (m  s-1 ) 

V  local  potential  (V) 

Vcep  output  voltage  (V) 

W  cell  stack  width;  width  (mm) 

y2-  molar  fraction  of  species  i 

Yj  mass  fraction  of  species  i 

Greaks  symbols 
a[m  molar  coefficient 

P  transfer  coefficient 

8a  anode  thickness  (mm) 

8C  cathode  thickness  (pun) 

8e  electrolyte  thickness  (pun) 

A/ihyd-oxy  enthalpy  change  of  H2  oxidation  reaction 
(kj  mol-1 ) 

AGhyd-oxy  Gibbs  free  enthalpy  change  of  H2  oxidation  reac¬ 
tion  (kj  mol-1 ) 
s  porosity 

£sofc  fuel  efficiency 


r 

degree  of  flow  uniformity 

ii 

overpotential  (V) 

/Ceff 

effective  thermal  conductivity  (W  m-1 1<-1 ) 

Kf 

thermal  conductivity  of  multi-component  mixture 
(Wm-1 1<-1) 

K'i 

thermal  conductivity  of  species  i  (W  m-1  K-1 ) 

Ks 

thermal  conductivity  of  solid  material  (W  m-1  K-1 ) 

^air 

air  ratio 

Meff 

effective  dynamic  viscosity  (kgm-1  s-1 ) 

fii 

dynamic  viscosity  of  species  i  (kgm-1  s-1 ) 

V 

kinematic  viscosity  (m2  s-1 ) 

P 

density  (kgm-3) 

X 

tortuosity 

a 

heat  transfer  coefficient 

Superscripts 

in 

feed  conditions  (fuel  and  air  channel  inlet) 

Subscripts 

a 

anode 

act 

activation  overpotential  (V) 

air 

air  gas  channel;  air  feed  header 

ave 

average  value 

bulk 

bulk  average  in  the  cross-sectional  area  of  the  cell 
stack 

c 

cathode 

con 

concentration  overpotential  (V) 

C02 

carbon  dioxide 

e 

electrolyte 

f 

fluid 

fuel 

fuel  gas  channel;  fuel  feed  header 

H2;  hyd 

hydrogen 

H20 

steam 

hyd-oxy 

reaction  of  hydrogen  and  oxygen 

int 

interconnect 

k 

x,  y,  and  z  directions  in  Cartesian  coordinates 

n 

total  number  of  rib-channels 

02;  oxy 

oxygen 

rib 

rib-channel 

s 

solid 

TPB 

three-phase  boundary 

This  paper  is  to  develop  a  comprehensive  mathematical  model 
on  the  basis  of  our  original  studies  [27-33]  to  accomplish  the  sim¬ 
ulation  of  the  one-cell  stack  proposed  by  Yakabe  et  al.  [34].  The 
present  mathematical  model  cannot  only  be  used  to  describe  the 
physical-electrochemical  phenomena  in  the  unit  cell,  but  also  for 
the  one-cell  stack.  Moreover,  self-design  methods  are  applied  to 
solve  the  model.  Subroutines  programmed  by  VC++  are  used  to 
solve  the  electrochemical  model  because  only  the  thermo-fluid 
model  can  be  solved  in  the  commercial  software  FLUENT.  Finally, 
influences  of  the  non-uniformity  of  gas  flow  rates  on  the  thermal 
and  electrical  performance  of  the  one-cell  stack  are  analyzed. 

2.  Model  development 

2.1.  Model  description 

The  one-cell  stack  consists  an  interconnect  structure  and  a 
three-layer  region  composed  of  two  ceramic  electrodes,  anode  and 
cathode,  separated  by  a  dense  ceramic  electrolyte  (often  referred 
to  as  PEN  (Positive-electrode/Electrolyte/Negative-electrode)).  The 
typical  geometry  of  a  planar  SOFC  proposed  by  Yakabe  et  al.  [34] 
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Fuel  exhaust  header 


Air  exhaust/feed  header  for 
co- flow/counter- flow  condition 


(b)  Air  feed/exhaust  header  for 
co-flow/counter-flow  condition 


Fuel  feed  header 


Fig.  1.  The  SOFC  stack  (a)  proposed  by  Yakabe  [34]  and  the  one-cell  stack  in  this  study  (b). 


is  depicted  in  Fig.  la.  The  one-cell  stack  considered  in  this  study 
is  shown  in  Fig.  lb.  The  computational  model  is  composed  of 
nine  components:  (1)  fuel  feed  header,  (2)  fuel  exhaust  header,  (3) 
fuel-side  interconnect  with  12  rib-channels,  (4)  12  fuel  channels 
(FH1 -FH12,  the  fuel  channel  (FH)  on  the  right  side  is  defined  as  FH1 ), 
(5)  PEN  structure,  (6)  12  air  channels  (AH1-AFI12,  the  air  channel 
(AH)  on  the  right  side  is  defined  as  AH1),  (7)  air-side  interconnect 
with  12  rib-channels,  (8)  air  feed  header  and  (9)  air  exhaust  header. 
The  schematic  of  the  one-cell  stack  and  its  repeating  cell  unit  pre¬ 
sented  in  Fig.  lb  is  shown  in  Fig.  2.  In  this  model,  the  thicknesses 
of  anode,  cathode,  electrolyte  and  interconnect  are  0.5  or  1.0,  0.25, 
0.05  and  2.0  mm,  respectively. 

2.2.  Thermo-fluid  model 

The  physical-chemical  transport  phenomena  in  a  SOFC  sys¬ 
tem  are  strongly  coupled.  Hence,  these  phenomena  are  classified 
into  the  following  categories:  (1 )  mass  transfer  in  gas  feed/exhaust 
headers,  gas  channels  and  porous  electrodes;  (2)  heat  transfer  in 
all  constituent  materials;  (3)  electrochemical  reactions  in  the  lay¬ 
ers  next  to  the  interfaces  between  electrolyte  and  electrodes.  In 
this  study,  the  solid  and  fluid  domains  are  all  divided  into  some 
discrete  meshes.  As  for  each  computational  mesh,  the  conservation 
equations  of  mass  continuity,  momentum,  energy  and  species  are 
solved  using  the  commercial  software  FLUENT  code  based  on  the 
finite  volume  method. 


The  mass  continuity  equation  is  written  as 

V.(spV)  =  Sm  (1) 

where  s  and  V  are  the  porosity  of  porous  electrodes  and  the  velocity 
vector,  respectively.  The  governing  equation  is  not  only  applicable 
to  porous  electrodes,  but  also  gas  channels  including  fuel  and  air 
channels.  And  £  is  used  to  distinguish  the  porous  electrodes  and  gas 
channels,  in  the  channels  s  =  1.  Both  air  and  fuel  flows  are  considered 
as  ideal  gas  mixtures,  therefore,  the  effective  density  of  the  multi¬ 
components  gas  mixture  based  on  ideal  gas  law  (IGL)  is  given  as 
follows  [31]: 


where  p,  R  and  T  are  the  pressure,  the  universal  gas  constant  and 
temperature,  respectively.  Yk  and  Mk  are  the  mass  fraction  and 
molecular  weight  of  species  k,  respectively.  Sm  is  the  source  term, 
which  describes  the  change  of  species  mass  caused  by  electro¬ 
chemical  reactions.  In  this  paper,  the  electrochemical  reactions  is 
assumed  that  only  occur  in  the  interfaces  between  the  electrolyte 
and  electrodes,  so  in  the  fuel  and  air  channels,  £  =  1  and  Sm  =  0.  And 
at  the  reaction  layers  next  to  the  interfaces  between  electrolyte  and 


Interconnect 


■  4mm- 

x 

(b) 


Fig.  2.  Schematic  of  the  one-cell  stack  (a)  and  the  cross-sectional  structure  of  the  cell  unit  (b). 
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Table  1 

Coefficients  of  Sutherland’s  viscosity  law  for  different  gas  components  [35]. 


Gas  species 

fio  (kgm-1  s-1) 

To  (K) 

S(I<) 

h2 

8.411E-6 

273.11 

96.67 

h2o 

1.703E-5 

416.67 

861.11 

co2 

1.370E-5 

273.11 

222.22 

o2 

1.919E-5 

273.11 

138.90 

n2 

1.663E-5 

273.11 

106.67 

porous  electrodes  Sm  is  as  follows: 

Sm  =  Yji  i  £  |H2,  H20  and  02}  (3) 

i 

where  is  the  mass  flux  of  species  i  during  the  electrochemical 
reaction  process. 

The  momentum  equation  is  as  follows: 


V  •  (epW)  =  -sVp  +  V  •  (e/XeffW)  +  Sv 


(4) 


where  /xeff  is  the  effective  dynamic  viscosity  of  the  mixture  gas,  and 
is  calculated  by  using  ideal  gas  mixing  law  based  on  kinetic  theory 
(MixKin)  [30]: 


A^eff  —  ^  ^ 


y^i 


f  'Yjifak 


(5) 


0  =  [1  +( Mt/^fc)1/2(Mfc/Mt)1/4]2  (g) 

*  [8(1  +  (Mj/Mj,))]1^2 

where  yz  is  the  mole  fraction  of  species  i.  Mz  and  Mk  are  molecular 
weight  of  species  i  and  k,  respectively.  /xz  and  /iik  are  dynamic  vis¬ 
cosities  of  species  i  and  k ,  respectively.  In  this  study,  the  dynamic 
viscosities  of  gas  components  are  defined  as  a  function  of  temper¬ 
ature  by  using  Sutherland’s  viscosity  law,  because  the  simulation 
of  the  SOFC  stack  involves  heat  transfer.  Sutherland’s  viscosity  law 
results  from  the  kinetic  theory  by  Sutherland  (1893)  using  an  ide¬ 
alized  intermolecular-force  potential.  Sutherland’s  law  with  three 
coefficients  is  [35]: 

/  r \3/2t0+s  ... 

M  >L°  \  To  j  T  +  S  (  ^ 

where  /x0  is  the  reference  value  of  viscosity.  T0  is  a  reference  tem¬ 
perature,  and  S  is  an  effective  temperature,  the  Sutherland  constant, 
which  is  characteristic  of  the  gas.  The  parameters  of  Sutherland’s 
law  for  gas  components  are  listed  in  Table  1. 

The  momentum  source  in  the  porous  electrodes,  Sv,  is  calculated 
by  the  following  formula: 

Sv  =  -(nr^+£/0BVk  M)  (8) 

where  I<  is  the  porous  electrode  permeability,  and  vk  is  the  velocity 
component  in  the  k  direction.  The  first  term  on  the  right  side  of 
Eq.  (8)  accounts  for  the  linear  relationship  between  the  pressure 
gradient  and  flow  rate  on  the  basis  of  the  Darcy’s  law.  The  second 


term  is  the  Forchheimer  term  concerning  the  inertial  force  effects 
[36].  In  the  fuel  and  air  channels,  s  =  1  and  Sv  =  0. 

In  general,  gas  species  transfer  mainly  by  convection  in  the  flow 
channels  and  diffusion  in  the  porous  electrodes.  The  species  con¬ 
servation  equation  is 


V  •  (spYiV)  =  V  •  (pD^effVYf)  +  Sz 


(9) 


where  Yz  is  the  mass  fraction  of  species  i.  Dz>eff  is  the  effective  gas 
diffusion  coefficient  of  species  i  in  the  porous  electrodes  and  fluid 
channels.  Based  on  the  dusty-gas  model  (DGM)  [34,37],  Dzeff  in  the 
porous  electrodes  can  be  expressed  as  follows: 


f  1  Qqm Vi  .  _J_A 

y  ^i,  m  ^Ki  J 


(10) 


where  r  is  the  tortuosity,  Di  m>  yz  and  DI<Z  are  the  molecular  diffusion 
coefficient,  the  molar  fraction  and  the  Knudsen  diffusion  coefficient 
of  the  species  i,  respectively.  Here,  a[m  is  defined  as  follows: 


where  Mz  is  the  molecular  weight  of  species  i,  and  Mm  is  average 
molecular  weight.  Knudsen  diffusion  occurs  in  the  porous  layer 
with  small  pores  or  under  low  pressure  when  the  mean  free-path 
of  molecules  is  larger  than  the  pore  size,  and  the  molecules  collide 
with  the  walls  more  often  than  between  themselves.  The  Knud¬ 
sen  diffusion  coefficient  for  component  i  in  the  multi-components 
mixture  is  calculated  based  on  the  free  molecule  flow  theory 


(12) 


where  dp  is  the  average  pore  size  of  the  electrodes.  In  a  multi¬ 
components  gas  system,  the  molecular  diffusion  coefficient  of  the 
component  i  is  given  by 


A',m  — 


i  -yt 

^  jCVk/^ik) 


(13) 


where  Dik  is  the  binary  diffusion  coefficient  in  the  system  with  gas 
components  i  and  k.  Sz  is  the  mass  production/consumption  rate  of 
species  i  depending  on  electrochemical  reactions.  The  production 
and/or  consumption  of  reactants  in  a  fuel  cell  are  proportional  to 
the  electronic  current  produced  by  the  electrochemical  reaction. 
For  hydrogen,  oxygen  and  water  species,  the  source  and/or  sink 
term  in  Eq.  (9)  can  be  expressed  as  follows: 


Sh2  =  ~2p  Mh2,  So2  =  -  — Mq2,  Sh2o  =  2p ^h2o  (14) 

where  i  is  the  charge-transfer  current  density.  F  is  the  Faraday  con¬ 
stant.  In  the  fuel  and  air  channels,  e  =  1  and  5Z  =  0. 

The  energy  equation  can  be  expressed  as 


V  •  (pCp VT)  =  V  •  (/ceffVT)  +  Sj  (15) 

where  Cp,  /ceff  and  Sj  are  the  specific  heat  capacity  of  composition- 
dependent  gas  mixture,  the  effective  thermal  conductivity  and  the 
energy  source  term,  respectively.  The  specific  heat  capacity  of  gas 


Table  2 

Polynomial  coefficients  for  specific  heat  capacity  of  gas  species  [35]. 


H20 

h2 

co2 

n2 

o2 

a0 

1.93780E+3 

1.4147E+4 

5.35446E+2 

1.02705E+3 

8.76317E+2 

a\ 

-1.18077 

1.7372E— 1 

1.27867 

2.16182E-2 

1.22828E-1 

a2 

3.64357E-3 

6.9E-4 

— 5.46776E— 4 

1.48638E— 4 

5.58304E-4 

a3 

-2.86327E-6 

- 

-2.38224E-7 

-4.48421  E-8 

-1.20247E-6 

a4 

7.59578E-10 

- 

1.89204E-10 

- 

1. 14741 E-9 

a5 

- 

- 

- 

- 

-5.12377E-13 

06 

- 

- 

- 

- 

8.56597E-17 
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Table  3 

Property  parameters  for  different  components  of  the  cell  unit  [6,30-34,37,39]. 


mixture  is  calculated  as  an  average  mass  fraction  of  the  pure  species 
based  on  the  mixing  law  [35]: 

Cp  =  ^>iCp,i  (16) 

i 

where  Cp  l  is  the  specific  heat  capacity  of  the  individual  species  i 
in  the  gas  mixture.  The  specific  heat  capacity  of  each  species  is 
calculated  as  a  function  of  temperature: 

m 

cP,i = Y,akJk  (17) 

k= 0 

where  the  polynomial  coefficients,  ak ,  proposed  by  Rose  and  Cooper 
[35]  are  listed  in  Table  2.  Heat  transfer  between  the  fluid  and  solid 
materials  is  limited  to  conduction  and  convection.  The  effect  of  radi¬ 
ation  is  neglected  in  this  study  because  it  is  very  small  relative  to 
the  other  kinds  of  heat  transfer.  Additionally,  the  effective  thermal 
conductivities  of  porous  electrodes  are  calculated  by  the  following 
equation  [31,38]: 


/Ceff  =  £/Cf  +  (l  -S)Ks 


(18) 


where  /Cf  and  /cs  are  thermal  conductivities  of  fluid  and  solid, 
respectively.  The  composition-dependent  thermal  conductivity  for 
multi-components  mixture  is  calculated  by  using  ideal  gas  mix¬ 
ing  law  based  on  kinetic  theory: 


Wi 

EiMik 


(19) 


where  (/)ik  is  defined  in  Eq.  (6).  K(  is  thermal  conductivity  of  species 
i.  The  thermal  conductivity  for  each  species  is  defined  using  kinetic 
theory  as 


15  R 

K;  =  —  —  l±i 

1  4 


4  CpjMi  1 
15  R  +  3 


(20) 


The  energy  source  term,  ST,  mainly  consists  of  reactions  to 
release  heat  and  ohmic  heat.  Based  on  the  local  current  density, 
thermodynamic  heat  on  the  anode  (except  the  ohmic  heat)  gener¬ 
ated  in  the  hydrogen  oxidation  reaction  is  calculated  as  follows: 


^T-thermo-a 


=  l 


Ah 


hyd-oxy 


2  F 


(^hyd-oxy  ^act  Peon) 


(21) 


where  the  subscripts  “hyd-oxy”  indicates  the  hydrogen  oxidation 
reaction.  Ahhyd_oxy  and  Ehyd-oxy  are  the  enthalpy  change  and  elec¬ 
tromotive  force  (EMF)  for  hydrogen  oxidation  reaction,  respectively. 
i) art  and  r] con  are  activation  and  concentration  overpotentials  of  the 
electrodes,  respectively.  Ohmic  heat  in  the  electrolyte  or  in  the  elec¬ 
trodes  is  also  locally  calculated  based  on  the  ohmic  resistivities  and 
current  densities.  The  ohmic  heat  sources  on  the  anode,  cathode 
and  electrolyte  are  calculated  as  follows: 

^T-ohmic-a  =  i^Rj-aSj_a 

<  ^T-ohmic-e  —  ^Rj-e^j-e  (22) 

k  ^T-ohmic-c  —  ^^j-c^j-c 

where  Rj_a,  Rj.ei  and  Rj.c  are  ohmic  resistivities  of  anode,  elec¬ 
trolyte  and  cathode,  respectively.  8j.a,  8j.a  and  8j_a  are  the  thickness 


of  anode,  electrolyte  and  cathode,  respectively.  Hence,  the  heat 
sources  on  the  anode,  electrolyte  and  cathode  are  as  follows: 


^T-a  =  i 


Ah 


hyd-oxy 


2F 


—  (^hyd-oxy  —  hact  ~  Peon) 


+  i2Rj.a8j_, 


ST-e  =  i2Rj-e8j_e 
ST-C  =  i2Rj-c8j-c 


(23) 


2.3.  Materials  and  property  settings 

The  anode  (Ni/YSZ  cermet)  and  cathode  (doped  lanthanum  man- 
ganite,  LaMn03)  are  modeled  as  porous  media.  The  electrolyte 
(yttria-stabilized  zirconia,  Zr02-8mol%  Y2O3)  and  interconnect 
(magnesium-doped  lanthanum  chromium  oxide  (LaMgCrOs))  are 
modeled  as  solid.  The  interface  between  electrodes  and  electrolyte 
are  treated  as  the  three-phase  boundary  (TPB).  Firstly,  all  the  faces 
in  electrodes  and  electrolyte  are  numbered.  If  the  coordinates  of 
face  center  is  same,  the  two  faces  are  mapped  and  made  as  the 
interface  between  electrodes  and  electrolyte.  Secondly,  the  each 
species  of  gas  are  numbered  and  to  obtain  the  specious  distribu¬ 
tions  in  the  interface.  Finally,  the  thermo-fluid  and  electrochemical 
models  are  solved.  In  this  paper,  solid  material  properties  used  in 
this  simulation  are  given  in  Table  3.  The  transport  properties  for 
the  computational  control  volume  in  the  present  study  are  listed  in 
Table  4. 


2.4.  Boundary  conditions 

The  inlet  flow  rate  into  fuel  and  air  feed  headers,  the  inlet  tem¬ 
perature  and  pressure,  species  concentration  at  fuel  and  air  feed 
inlet,  and  output  voltage  are  all  summarized  in  Table  5.  At  the 
exhaust  outlets  of  fuel  and  air,  the  fixed  pressure  boundary  condi¬ 
tion  is  adopted  and  fixed  as  1  atm.  The  thermal  boundary  conditions 
at  the  top  and  bottom  surfaces  of  the  one-cell  stack  are  regarded 
as  thermally  adiabatic.  As  it  is  difficult  to  determine  the  exact 
boundary  condition  at  the  other  edge  parts,  the  adiabatic  bound¬ 
ary  conditions  are  employed  at  the  left  and  right  edges.  For  the 
boundary  between  the  solid  and  the  fluid,  the  following  continuity 
condition  is  imposed  [34]: 

[/ceff(x)VTs(x)]  x  n  =  a[T{(x)  -  Ts(x)]  (24) 

where  /ceff(x)  is  the  thermal  conductivity  at  x  on  the  boundary,  a 
is  the  heat  transfer  coefficient,  n  is  the  unit  vector  normal  to  the 
boundary,  Ts  and  Tf  are  temperatures  of  the  solid  and  fluid  at  x  on 
the  boundary,  respectively. 


Table  4 

Transport  property  setting  for  computational  control  volume  in  the  present 
simulation. 


p  (kg  ITT3) 

Ke  fKWm-MC1) 

CpCJkg-’K-’) 

/^eff  (kgm-1  s-1) 

D,-eff  (m2  s-1) 

IGL 

Eq*  (18) 

Mixing  law 

MixKin 

DGM 

IGL:  ideal  gas  law;  MixKin:  ideal  gas  mixing  law  based  on  kinetic  theory;  DGM: 
dusty-gas  model. 
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Table  5 

Parameters  and  conditions  in  this  study  [6,23-26,39]. 


Description 

Value 

Species  molar 

Hydrogen,  yH2 

8/11 

concentration  at 

Carbon  dioxide,  yCo2 

2/11 

fuel  feed  inlet 

Water,  yH2o 

1/11 

Total 

1.0 

Species  molar 

Oxygen,  yo2 

0.21 

concentration  at 

Nitrogen,  yN2 

0.79 

air  feed  inlet 

Total 

1.0 

Exchange  current 

Anode,  i0i a 

5300 

(Am-2) 

Cathode,  i0>c 

2000 

Inlet  delivery  rate 

Fuel  feed  header,  Vfuei 

0.6/1.8 

(ms-1) 

Air  feed  header,  Vair 

4.0/6.0 

Inlet  temperature 

Fuel  feed  header,  Tfuel 

625/675 

(°C) 

Air  feed  header,  Tair 

625/675 

Operating  voltage  (V) 

Vcdl 

0.6 

Operating  pressure  (kPa) 

P 

101.3 

Porosity  (%) 

Anode,  ea 

0.45/0.25 

Cathode,  ec 

0.45/0.25 

Tortuosity 

Anode,  ra 

3.0 

Cathode,  rc 

2.0 

Pore  size  of  cathode 

Anode,  dp,a 

1.0/3.0 

(pan) 

Cathode,  dp>c 

3.0 

Length  (mm) 

Anode,  Za 

50 

Cathode,  lc 

50 

Electrolyte,  le 

50 

Thickness  (pm) 

Cathode,  8C 

250 

Electrolyte,  8e 

50 

Anode,  <5a 

500/1000 

Cell  unit  number  of  the  stack 

n 

12 

Width  (mm) 

Stack,  W 

46 

Fuel  channel  in  the  cell  unit,  Wm 

2.0 

Air  channel  in  the  cell  unit,  WA H 

2.0 

Rib-channel,  Wrib 

2.0 

Height  (mm) 

Fuel  channel,  dFH 

1.0 

Air  channel,  dAH 

1.0 

Rib-channel,  drib 

1.0 

Diameter  of  gas  inlet  (mm) 

Fuel  feed/exhaust  header  and  air 
feed/exhaust  header 

4.0 

2.5.  Electrochemical  model 


During  the  energy  transforming  process,  when  the  charge  trans¬ 
fer  reaction  at  the  electrolyte-electrode  interface  is  too  slow  to 
provide  ions  at  the  rate  required  by  the  demand  of  current,  the 
activation  polarization  occurs  and  is  defined  by  the  Butler-Volmer 
equation  [36]: 


i  =  i*o  jexp  (  - 

-exp(l 

-/3)^»?act]} 

(29) 

Then, 

2RT  •  u-l  1 

^act.a  =  Sinh  1 

(2l0,a; 

) 

(30) 

^ct.c=  n^smh  1 

\  2l0,c 

) 

(31) 

where  /3  is  the  transfer  coefficient  (0  <  /3  <  1 ),  the  transfer  coefficient 
is  considered  to  be  the  fraction  of  change  in  polarization  that  leads 
to  a  change  in  reaction  rate  constant,  and  its  value  is  usually  0.5  for 
the  fuel  cell  application  [39],  so  p  «  0.5  in  this  simulation,  ne  is  the 
number  of  electrons  participating  in  the  reaction,  i0  is  the  exchange 
current  densities  at  the  electrodes.  i0,a  and  i0iC  are  the  exchange 
current  densities  at  the  anode  and  the  cathode,  respectively.  In  this 
study,  exchange  current  densities  at  the  anode  and  the  cathode  are 
5300  Am-2  and  2000  Am-2,  respectively. 

The  concentration  overpotential  is  defined  as  the  difference 
between  the  EMF  calculated  based  on  the  mean  concentration  over¬ 
potential  Ebuik  and  that  based  on  the  local  concentration  on  the  TPB 
Etpb 


_  _  RT 

Peon  =  Ebuik  _  Etpb  =  2f  11 


Po2,bulk  PH2,bulk  x  Ph20,TPB 

dJ/2™  Ph2, TPB  xPh20, bulk 

^02,TPB  z  2  ’ 


(32) 


where  Pn2,buik  and  PH2o,buik  are  the  hydrogen  and  water  partial 
pressures  in  fuel  channel,  respectively,  andpH2jpB  andpH2o,TPB  are 
those  at  the  interface  of  electrolyte/anode,  respectively.  p02,  bulk  and 
Po2  ,tpb  are  oxygen  partial  pressure  in  air  channel  and  at  the  interface 
of  electrolyte/cathode. 

The  local  current  density,  z,  and  the  ohmic  losses,  pohm,  of  the 
cell  components  are  calculated  by  using  Ohm’s  and  Kirchhoff  s  law 
with  the  values  of  ohmic  resistivity,  Rj ,  as  shown  in  Table  3: 


'Johm  =  ^2iRj  J  e  fa.  c,  e  and  int} 
i 


(33) 


2.5.1.  Reactions  description 

The  oxidant  reduction  reaction  occurring  at  the  cathode  is 
expressed  as  follows: 

j02  +  2e~  -»■  O2-  (25) 


The  details  of  individual  overpotentials  above  can  be  obtained 
from  Suzuki  et  al.  [38]  and  Chan  et  al.  [39]. 

Finally,  the  following  local  relations  exist  among  the  electromo¬ 
tive  force  F,  overpotential  loss  p,  current  density  z,  and  the  output 
voltage  of  the  cell  stack,  Vcen: 


The  oxygen  ions  transfer  through  the  electrolyte  and  then  into 
the  active  reaction  layers  of  anode.  The  electrochemical  reaction  of 
fuel  at  the  anode  is 


H2  +  O2-  -»■  H20  +  2e~ 

(26) 

Hence,  the  overall  reaction  is 

H2  +  102  -y  h2o 

(27) 

2.5.2.  Electrochemical  reactions  dynamics 

According  to  the  Faraday  law,  the  reaction  rates  depend  on  the 
current  density  z  [31  ]: 

d  f  dO? 

i  =  2F-j-  =  4F 

dt  dt 

(28) 

where  df/dt  and  d02  /dt  are  the  molar  consumption  rates  of  fuel  and 
oxygen  at  the  anode  and  the  cathode,  respectively. 


Ehyd-oxy  ( Pact  +  Peon  +  Pohm)  —  ^cell  (34) 
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Table  6 

Performance  factors  of  the  cell  stack. 


Fuel  utilization  factor 
(36)t/fuei  =  'TneLW 

2%n/fuel 

Air  ratio 

Power  density 

(38) PsoFC  =  lave^cell 

Fuel  efficiency 

(39) Ssofc  =  ‘7^"™ 

)'H2QH2"f“1 


where  iave  is  the  average  total  current 
density,  I  and  W  are  the  cell  length  and 
width,  nair  and  Hfuei  are  air  and  fuel  feed 
header  inlet  flow  rates,  Q.  is  lower  heat¬ 
ing  value.  The  superscript  in  denotes 
feed  conditions  (fuel  and  air  feed  header 
inlet). 


where  Ehyd-oxy  is  computed  by  Nernst  equation: 

1/2  \ 


AG 


^hyd-oxy 


rhyd-oxy 


2  F 


RT _ 

+  2 F  n 


Ph2Pq2 
Ph,o 


(35) 


where  AGhyd_oxy  is  the  Gibbs’s  free-energy  change  of  the  hydrogen 
oxidation  reaction,  Ph2o>  Po2  and  Ph2  are  the  partial  pressures  of 
water  gas,  oxygen  and  hydrogen,  respectively.  The  partial  pressure 
of  species  t,  pz,  is  obtained  as  a  product  of  the  total  local  pressure 
and  the  molar  fraction  yz. 

If  the  output  voltage  of  the  cell  stack  Vcell  is  given,  the  local  cur¬ 
rent  density  i  can  be  obtained  from  these  relations  above,  and  vice 
versa.  In  the  present  model,  the  cell  terminal  voltage  defined  as 
Vceii  is  given  as  a  signal  to  start  the  computation  for  different  oper¬ 
ating  conditions  (called  operating  voltage),  and  thereby  both  the 
current  densities  and  potential  fields  are  calculated  based  on  an 
equivalent  electrical  circuit  model  illustrated  in  Fig.  3.  In  the  one¬ 
cell  stack,  current  is  assumed  to  pass  perpendicular  to  the  gas  flow 
direction  through  the  anode,  electrolyte,  cathode  and  interconnec¬ 
tor  (including  ribs).  And  this  stack  is  divided  into  large  number  of 
the  meshes.  For  the  meshes  with  same  centric  coordinate  in  the 
PEN  and  interconnector  (including  ribs),  the  meshes  are  used  to 
the  circuit  equivalent.  Since  the  steady-state  performance  of  the  cell 
stack  is  discussed  in  the  present  study,  the  capacitance  between  the 
electrolyte  and  the  electrodes  is  neglected. 


2.6.  Model  performance  factors 


The  fuel  utilization  factor  is  the  fraction  of  the  total  inlet  fuel  con¬ 
sumed  to  produce  electricity  in  the  cell  stack.  The  air  ratio  reflects 
the  excess  air,  with  respect  to  the  stoichiometrically  demand,  oth¬ 
erwise  the  excessive  will  cooling  the  cell  stack.  Table  6  presents  the 
mathematical  definition  of  these  performance  factors. 


2.7.  Model  calculation  process 

The  governing  Eqs.  (1),  (4),  (9)  and  (15),  together  with  their 
relevant  boundary  conditions  and  initial  conditions,  are  solved 
with  computational  fluid  dynamics  software  FLUENT  by  finite  vol¬ 
ume  method  coupling  the  electrochemical  Eq.  (32).  The  modelling 
domain  consists  of  the  fuel  feed/exhaust  headers,  fuel  channel,  air 
feed/exhaust  headers,  air  channel,  interconnect  and  PEN  structure, 
as  shown  in  Fig.  1.  All  the  source  terms  in  the  mass  continuity, 
species  conservation  and  energy  equations  must  be  adopted  in 
the  electrochemical  model  because  this  model  is  coupled  with  the 
thermo-fluid  one.  Moreover,  the  electrochemical  model  must  be 
solved  by  the  subroutines  developed  in  this  study,  and  the  model 
is  not  only  originated  from  the  software  Fluent. 

The  finite-volume  Navier-Stokes  and  transport  equations  are 
solved  to  obtain  the  gas  species  concentrations,  velocities  and  tem¬ 
peratures  at  each  computational  control  volume  in  the  one-cell 
stack.  These  information  are  supplied  for  the  electrochemical  model 
to  calculate  the  local  current  density,  and  this  model  is  processed 
by  the  self-developed  subroutine  interface  using  the  user  defined 
function  (UDF)  provided  by  FLUENT.  Then,  the  resultant  current 
density  is  applied  to  obtain  the  hydrogen  reaction  rate,  heat  source 
and  species  sources.  Gas  species  concentrations,  velocities  and  tem¬ 
perature  distributions  are  calculated  for  the  next  iteration,  and  so 
on,  until  the  convergence  of  solution  is  achieved. 

All  the  cases  in  this  study  are  shown  in  Table  7,  and  the  rest 
operating  conditions  and  parameters  are  listed  in  Table  5. 


3.  Results  and  discussion 

3  A.  Effects  of  flow  patterns  on  the  cell  stack  performance  based 
on  flow  uniformity  analysis 

Fig.  4  shows  the  velocity  distributions  in  fuel  and  air  channels  for 
counter  flow  in  case  1  and  co-flow  in  case  2  (Table  7).  In  this  figure, 
U[  is  the  velocity  at  the  ith  rib-channel,  and  uaVe  is  the  averaged 
mean  velocity  of  all  the  channels.  The  flow  rates  at  the  channels 
close  to  the  entrance  and  exist  of  the  gas  feed/exhaust  headers  are 
higher  than  those  in  the  middle  channels  both  in  counter  flow  and 
co-flow  conditions. 

For  quantitative  estimation,  the  degree  of  flow  uniformity  is 

defined  as  follows  [6]:  r  =  1  -  {(l/n)^"=1((uf  -  uave)/uave)2}^2 , 
where  n  =  12  denotes  the  total  number  of  rectangular  rib-channels 
in  the  one-cell  stack. 

As  for  the  two  flow  configurations  calculated,  the  cell  stack 
operating  in  counter  flow  case  has  a  more  uniform  flow  veloc¬ 
ity  distribution  in  both  fuel  and  air  channels.  As  for  the  counter 
flow  condition  in  case  1,  the  degrees  of  flow  uniformity  in  fuel 
and  air  channels  are  0.9620  and  0.8985,  respectively.  Accordingly, 
the  degrees  in  fuel  and  air  channels  are  0.9289  and  0.8970  under 
co-flow  condition  in  case  2,  respectively. 


Table  7 

Computational  parameters  for  all  the  calculated  cases. 


Case 

1 

2 

3 

4 

5 

6 

7 

8 

9 

Flow  pattern 

Counter 

Co- flow 

Counter 

Counter 

Counter 

Counter 

Counter 

Counter 

Counter 

flow 

flow 

low 

flow 

flow 

flow 

flow 

flow 

Air  inlet  velocity,  Vair  (m  s-1 ) 

6.0 

6.0 

4.0 

6.0 

6.0 

6.0 

6.0 

6.0 

6.0 

Fuel  inlet  velocity,  Vfuei  (ms-1) 

0.6 

0.6 

0.6 

1.8 

1.8 

1.8 

1.8 

1.8 

1.8 

Air  inlet  temperature,  Tair  (°C) 

625 

625 

625 

625 

675 

625 

625 

625 

625 

Fuel  inlet  temperature,  Tfuei  (°C) 

625 

625 

625 

625 

625 

675 

625 

625 

625 

Porosities  of  anode  and  cathode,  e  (%) 

45 

45 

45 

45 

45 

45 

25 

45 

45 

Pore  size  of  the  anode,  dp,a  (fxm) 

3.0 

3.0 

3.0 

3.0 

3.0 

3.0 

3.0 

1.0 

3.0 

Anode  thickness,  <5a  (p,m) 

500 

500 

500 

500 

500 

500 

500 

500 

1000 

Values  in  bold  of  Cases  2-4  are  differences  from  those  of  the  case  1,  and  values  in  bold  of  cases  5-9  are  differences  from  those  of  the  case  4. 
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Fig.  4.  Velocity  distributions  in  fuel  (a)  and  air  channels  (b)  for  counter  flow  in  case  1  and  co-flow  in  case  2  listed  in  Table  7. 


Table  8 

Cell  stack  performance. 


Current  density  distributions  for  counter  flow  in  case  1  and 
co-flow  in  case  2  are  shown  in  Fig.  5,  and  their  temperature  dis¬ 
tributions  of  PEN  structure  are  shown  in  Fig.  6.  The  average  values 
of  voltage,  current  density  and  the  temperature  of  PEN  structure 
for  all  the  calculated  cases  are  itemize  in  Table  8  as  well  as  power 
density,  fuel  utilization  factor,  air  ratio  and  fuel  efficiency. 

In  the  counter  flow  case,  the  current  density  decreases  along 
the  fuel  flow  direction,  and  is  high  near  the  fuel  inlet  due  to  the 
high  temperature  leading  to  small  ohmic  overpotential  and  Nernst 
potential  in  that  region.  This  is  because  EMF  gradually  decreases 
from  fuel  inlet  toward  outlet  as  02  and  H2  are  consumed  and  H20 
is  produced.  In  contrast,  the  maximum  current  density  appears  in 
the  middle  of  the  cell  stack  in  the  co-flow  case.  This  is  because 
the  overpotentials  are  very  temperature  sensitive  and  there  is  a 
trade-off  between  the  temperature  increase  and  the  reactant  con¬ 


centration  decrease.  The  characteristics  of  current  density  in  both 
the  two  flow  cases  are  similar  to  those  of  the  planar  one-cell  stack 
with  uniform  flow  rates  in  fuel  and  air  inlets  [11].  In  addition,  the 
characteristic  in  co-flow  case  is  similar  to  that  of  the  planar  cell  unit 
[34,36].  Compared  with  the  co-flow  case,  the  current  density  distri¬ 
bution  in  counter  flow  case  is  more  uniform  in  the  width  direction 
of  the  cell  stack  due  to  the  better  uniform  fuel  and  air  flows.  The 
high  current  densities  appear  in  the  location  where  are  channels 
with  large  fuel  delivery  rate.  From  Table  8,  it  is  suggested  that  the 
average  current  density,  power  density,  fuel  utilization  factor  and 
fuel  efficiency  for  counter  flow  case  in  case  1  are  higher  than  those 
for  co-flow  case  in  case  2. 

The  average  temperature  of  PEN  structure  is  1057.2  K  in  counter 
flow  case  shown  in  Table  8.  Furthermore,  the  PEN  structure  tem¬ 
perature  rises  rapidly  in  the  air  flow  direction,  reaching  a  maximum 


Fig.  5.  Current  density  distributions  for  counter  flow  in  case  1  (a)  and  co-flow  in 
case  2  (b)  listed  in  Table  7. 


Fig.  6.  Temperature  distributions  of  PEN  structure  for  counter  flow  in  case  1  (a)  and 
co-flow  in  case  2  (b)  listed  in  Table  7. 
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exhaust 


Fig.  7.  Averaged  temperature  difference  of  PEN  structure  for  counter  flow  in  case  1 
and  co-flow  in  case  2. 


near  the  fuel  inlet,  and  then  gradually  drops.  Air  flow  has  most  effec¬ 
tively  cooling  effect  near  the  air  inlet,  and  carries  heat  generated  in 
the  PEN  structure,  and  thus  moves  the  hotspot  toward  the  air  out¬ 
let.  The  average  temperature  of  PEN  structure  is  1041.4  K  in  co-flow 
case.  However,  it  rises  uniformly  along  the  direction  of  fuel  flow, 
and  is  highest  near  the  fuel  outlet  close  to  fuel  exhaust  header,  as 
shown  in  Fig.  7.  The  temperature  gradient  of  PEN  structure  in  co¬ 
flow  case  is  smaller  than  that  in  counter  flow  case  from  air  inlet  to 
outlet.  This  is  due  to  the  offsetting  effects  of  air  near  the  inlet,  at  its 
coolest,  being  aligned  with  the  fuel  inlet. 

Compared  with  co-flow  case,  it  is  should  be  noted  that  the  tem¬ 
perature  distribution  of  PEN  structure  for  counter  flow  case  is  more 
uniform  in  the  width  direction  of  the  cell  stack.  This  is  because  the 
better  uniform  fuel  and  air  flows  leading  to  better  current  density 
distribution.  The  high  temperature  of  PEN  structure  is  intensively 
concentrated  in  the  regions  close  to  the  fuel  exhaust  header.  The 
temperature  characteristics  of  PEN  structure  in  counter  flow  case 
are  similar  to  those  of  planar  one-cell  stack  with  uniform  flow  rates 
in  fuel  and  air  inlets  [11  ].  In  addition,  the  characteristics  in  co-flow 
case  is  similar  to  the  relevant  findings  in  the  literature  [31,31,34,36]. 
The  results  in  counter  flow  case  are  also  similar  to  that  of  the  planar 
cell  unit  [30,31]. 

As  a  conclusion,  the  cell  stack  in  counter  case  has  a  better  per¬ 
formance  (such  as  average  current  density,  power  density,  fuel 
utilization  factor,  fuel  efficiency  and  so  on)  than  that  in  co-flow  case. 
To  achieve  a  more  uniform  PEN  structure  temperature  distribution, 
a  deeply  investigation  into  the  effect  of  operating  conditions,  phys¬ 
ical  and  structural  parameters  of  the  cell  stack  on  its  performance 
based  on  flow  uniformity  analysis  in  counter  case  is  described  in 
the  follows. 
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Fig.  9.  Temperature  distributions  of  PEN  structure  for  counter  flow  in  case  1  (a),  case 
3  (b)  and  case  4  (c)  with  different  air  inlet  flow  rates. 


3.2.  Effects  of  inlet  velocity  on  the  cell  stack  performance  based 
on  flow  uniformity  analysis 

The  velocity  distributions  in  fuel  and  air  channels  and  tempera¬ 
ture  distributions  of  PEN  structure  along  the  air  flow  direction  for 
cases  1, 3  and  4  are  illustrated  in  Figs.  8-10,  as  well  as  the  details  of 
average  temperature  difference  of  PEN  structure. 

As  shown  in  Fig.  8,  the  degrees  of  flow  uniformity  in  different 
fuel  channels  decrease  from  0.9620  to  0.9537  and  0.9429  with  the 
decrease  of  air  inlet  velocity  and  the  increase  of  fuel  inlet  velocity. 
Meanwhile,  the  degrees  of  flow  uniformity  in  air  channels  slightly 
increase  from  0.8985  to  0.9034  and  0.9010  when  the  air  inlet  veloc¬ 
ity  decreases  and  the  fuel  inlet  velocity  increases.  Consequently, 
fuel  flows  are  less  uniform  with  the  decrease  of  air  inlet  velocity 
and  the  increase  of  the  fuel  inlet  velocity. 
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Fig.  8.  Velocity  distributions  in  fuel  (a)  and  air  channels  (b)  for  counter  flow  in  cases  1,  3  and  4  listed  in  Table  7. 
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Fig.  10.  Averaged  temperature  difference  of  PEN  structure  for  counter  flow  in  cases 
1,  3  and  4  listed  in  Table  7. 


As  listed  in  Table  8,  the  average  temperature  of  PEN  struc¬ 
ture  rises  from  1057.2  K  to  1148.6 1<  when  the  air  inlet  velocity 
decreases  from  6.0ms_1  to  4.0  ms-1.  The  average  current  density 
increases  from  4156.2  Am-2  to  4267.4 Am-2  with  the  decrease  of 
air  inlet  velocity  due  to  the  ohmic  overpotential  decrease  with  the 
increase  of  PEN  structure  temperature.  Moreover,  the  power  den¬ 
sity,  fuel  utilization  factor  and  fuel  efficiency  all  increase  with  the 
decrease  of  air  inlet  velocity.  It  is  to  be  noted  that  the  decrease 
of  air  ratio  is  the  reason  why  the  average  temperature  of  PEN 
structure  rises  with  the  decrease  of  air  inlet  velocity.  The  aver¬ 
age  temperature  of  PEN  structure  rises  from  1057.2  K  to  1081.9  K 
when  the  fuel  inlet  velocity  increases  from  0.6  ms-1  to  1.8  ms-1, 
which  is  because  more  heat  produced  by  H2  oxidation  reaction. 
The  reaction  is  promoted  by  the  increase  of  average  current  den¬ 
sity  with  enhanced  electrochemical  reactions  by  the  more  fuel 
supply.  With  the  increase  of  fuel  inlet  velocity,  the  power  density 
increases,  but  the  air  ratio,  fuel  utilization  factor  and  fuel  efficiency 
decrease. 

It  can  be  observed  in  Fig.  9  that  the  distribution  of  PEN  structure 
temperature  is  less  uniform  with  the  decrease  of  air  inlet  velocity 
and  the  increase  of  the  fuel  inlet  velocity  due  to  the  worse  unifor¬ 
mity  of  fuel  flows.  The  temperature  difference  of  PEN  structure  in 
the  width  direction  of  the  cell  stack  increases  with  the  decrease  of 
air  inlet  velocity,  and  the  high  temperature  region  expands  with  the 
increase  of  fuel  inlet  velocity. 

Fig.  10  presents  the  average  temperature  difference  of  PEN 
structure  in  the  air  flow  direction  for  cases  1,  3  and  4.  The  tem¬ 
perature  gradient  of  PEN  structure  increases  when  air  inlet  velocity 
decreases  or  fuel  inlet  velocity  increases. 


Fig.  12.  Temperature  distributions  of  PEN  structure  for  counter  flow  in  case  4  (a), 
case  5  (b)  with  higher  air  inlet  temperature  and  case  6  (c)  with  higher  fuel  inlet 
temperature  listed  in  Table  7. 

3.3.  Effects  of  inlet  temperature  on  the  cell  stack  performance 
based  on  flow  uniformity  analysis 

Velocity  distributions  in  fuel  and  air  channels  and  temperature 
distributions  of  PEN  structure  are  illustrated  in  Figs.  12  and  13,  as 
well  as  the  details  about  the  average  temperature  difference  of  PEN 
structure  along  the  air  flow  direction  for  cases  4-6. 

As  shown  in  Fig.  11,  the  degree  of  flow  uniformity  decreases 
from  0.9429  to  0.9316  in  fuel  channels,  and  increases  from  0.9010 
to  0.9104  in  air  channels  with  the  increase  of  air  inlet  tempera¬ 
ture.  On  the  contrary,  with  the  increase  of  fuel  inlet  temperature, 
the  degree  in  fuel  channels  increases  from  0.9429  to  0.9537,  and 
decreases  from  0.9010  to  0.8967  in  air  channels.  Therefore,  air  flows 
are  more  uniform  with  the  increase  of  air  inlet  temperature.  Oth- 
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Fig.  11.  Velocity  distributions  in  fuel  (a)  and  air  channels  (b)  for  counter  flow  in  cases  4-6  listed  in  Table  7. 
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Fig.  13.  Averaged  temperature  difference  of  PEN  structure  for  counter  flow  in  cases 
4-6  listed  in  Table  7. 


erwise,  they  are  less  uniform  with  the  increase  of  the  fuel  inlet 
temperature. 

As  listed  in  Table  8,  the  average  temperatures  of  PEN  struc¬ 
ture  rise  from  1081.9  K  to  1137.3  K  and  1090.7  K,  when  the  air  and 
fuel  inlet  temperatures  increase  from  625  °C  to  675  °C,  respec¬ 
tively.  The  average  current  densities  increase  from  4873.8  Am-2 
to  5009.7 Am-2  and  4883.1  Am-2  with  the  increase  of  air  and 
fuel  inlet  temperatures,  respectively,  which  is  because  the  ohmic 
overpotential  decreases  due  to  the  increase  of  PEN  structure  tem¬ 
perature.  Moreover,  the  power  density,  fuel  utilization  factor  and 
fuel  efficiency  increase  with  the  increase  of  air  and  fuel  inlet  tem¬ 
peratures. 

Fig.  12  shows  that  the  distribution  of  PEN  structure  tempera¬ 
ture  is  more  uniform  in  the  width  direction  of  the  cell  stack  with 
the  increase  of  air  inlet  temperature.  This  case  is  due  to  the  better 
uniform  air  flows.  With  the  increase  of  fuel  inlet  temperature,  the 
temperature  difference  is  almost  constant  in  the  width  direction  of 
the  cell  stack. 

Fig.  13  presents  the  average  temperature  difference  of  PEN 
structure  in  the  air  flow  direction  for  cases  4-6  in  details.  The  tem¬ 
perature  gradient  of  PEN  structure  increases  with  the  individual 
increase  of  air  and  fuel  inlet  temperatures.  It  is  to  be  noted  that 
the  decrease  of  air  ratio  and  the  increase  of  average  current  density 
are  the  reasons  why  the  temperature  gradients  of  PEN  structure 
increase  with  the  increase  of  air  and  fuel  inlet  temperatures.  In 
the  higher  air  inlet  temperature  case  (case  5),  the  ascending  of  the 
average  temperature  gradient  is  2.49  K  mm-1,  and  the  descending 
is  1.77  Knurr1.  If  the  cell  stack  operates  on  the  condition  that  the 
ascending  of  average  temperature  gradient  is  less  than  its  descend¬ 
ing,  more  uniform  PEN  structure  temperature  distribution  may  be 
achieved  in  the  width  direction  of  the  cell  stack,  as  well  as  higher 
power  output,  fuel  utilization  factor  and  fuel  efficiency. 


T(K) 


Fig.  15.  Temperature  distributions  of  PEN  structure  for  counter  flow  in  case  4  (a), 
case  7  (b)  with  lower  porosities  of  electrodes,  case  8  (c)  with  smaller  pore  size  of 
anode  and  case  9  (d)  with  thicker  anode  thickness  listed  in  Table  7. 


3.4.  Effects  of  physical  and  structural  parameters  on  the  cell  stack 
performance  based  on  flow  uniformity  analysis 

Velocity  distributions  in  fuel  and  air  channels  and  tempera¬ 
ture  distributions  of  PEN  structure  along  the  air  flow  directions 
are  shown  in  Figs.  14-16,  as  well  as  the  details  about  the  average 
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Fig.  14.  Velocity  distributions  in  fuel  (a)  and  air  channels  (b)  for  counter  flow  in  cases  4,  7-9  listed  in  Table  7. 
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temperature  difference  of  PEN  structure  for  cases  4,  7-9,  respec¬ 
tively. 

As  shown  in  Fig.  14,  the  degree  of  flow  uniformity  in  fuel  chan¬ 
nels  decreases  from  0.9429  to  0.9392,  and  increases  from  0.9010 
to  0.9027  in  air  channels  with  the  decrease  of  anode  porosity.  The 
degrees  increase  from  0.9429  and  0.9010  to  0.9439  and  0.9019  with 
the  decrease  of  anode  pore  size  in  fuel  and  air  channels,  respectively. 
The  degree  increases  from  0.9429  to  0.9715  in  fuel  channels,  and 
decreases  from  0.9010  to  0.8995  in  air  channels  with  the  increase 
of  anode  thickness.  In  conclusion,  air  flows  are  more  uniform  with 
the  decrease  of  anode  porosity  and  its  pore  size,  but  are  less  uniform 
with  the  increase  of  anode  thickness. 

As  listed  in  Table  8,  the  average  temperatures  of  PEN  struc¬ 
ture  decline  from  1081.9  K  to  1076.4  K,  1078.1  K  and  1069.3  K 
when  the  anode  porosity  and  pore  size  decrease,  and  the  anode 
thickness  increases.  The  average  current  densities  decrease  from 
4873.8 Am-2  to  4737.1  Am-2  and  4754.5 Am-2  because  of  the 
larger  concentration  overpotentials  with  the  decrease  of  anode 
porosity  and  the  increase  of  anode  thickness.  Furthermore,  the 
power  density,  fuel  utilization  factor  and  fuel  efficiency  decrease 
with  the  decrease  of  anode  porosity  and  the  increase  of  anode 
thickness.  The  average  current  density  increase  from  4873.8  A  m-2 
to  5013.8  Am-2  owing  to  the  smaller  concentration  overpotentials 
when  anode  pore  size  decreases.  To  sum  up,  the  power  density,  fuel 
utilization  factor  and  fuel  efficiency  increase  when  anode  pore  size 
decreases. 

As  shown  in  Fig.  15,  the  temperature  distribution  of  PEN  struc¬ 
ture  is  more  uniform  in  the  width  direction  of  the  cell  stack  with 
the  decrease  of  anode  porosity  due  to  the  better  uniform  of  air 
flows.  The  temperature  difference  of  PEN  structure  decreases  in  the 
width  direction  of  the  cell  stack  when  anode  porosity  decreases,  but 
increases  with  the  increase  of  anode  thickness.  However,  the  tem¬ 
perature  difference  of  PEN  structure  is  almost  constant  with  the 
decrease  of  anode  pore  size. 

Fig.  16  shows  the  average  temperature  difference  of  PEN  struc¬ 
ture  in  the  air  flow  direction  for  cases  4,  7-9.  The  temperature 
gradient  of  PEN  structure  falls  down  when  anode  porosity  and 
anode  thickness  increase  because  of  the  increase  of  air  ratio, 
whereas  it  is  almost  unchanged  with  the  variation  of  anode  pore 
size,  as  mentioned  in  Ref.  [36].  This  is  because  the  heat  transfer 
in  the  porous  anode  is  dominated  by  the  streamwise  heat  conduc¬ 
tion,  which  is  approximately  proportional  to  the  effective  thermal 
conductivity  defined  in  Eq.  (18)  when  the  cross-sectional  area  is 
constant.  Therefore,  the  temperature  distribution  is  independent 
of  anode  pore  size.  Moreover,  the  PEN  structure  temperature  con¬ 
tinuously  rises  along  the  air  flow  directions  in  the  thicker  anode 
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Fig.  17.  Comparisons  between  the  calculated  flow  velocity  distributions  and  the 
experimental  (uave  is  the  averaged  velocity  of  12  channels). 


case  (case  9)  because  the  lower  heat  transfer  in  thicker  anode  will 
weaken  the  cooling  effect  of  fuel  flows. 

3.5.  Validation 

An  important  parameter  used  to  control  flow  rates  for  the  fluid 
flow  channel,  the  hydraulic  Reynolds  number,  Re  =  uaveDh/v,  is 
defined,  where  uave  is  the  averaged  velocity  from  12  rib-channels, 
v  is  the  kinematical  viscosity  of  the  fluid,  and  Dh  =  2  Wd/(W  +  d)  is 
the  hydraulic  diameter.  Note  that  W  and  d  are  the  width  and  height 
of  the  rib-channels. 

In  order  to  validate  the  numerical  model  developed  in  this  paper, 
the  numerical  flow  velocity  data  are  compared  with  the  exper¬ 
imental  reported  by  Huang  et  al.  [6].  The  comparison  between 
the  calculated  flow  velocity  distributions  and  the  experimental  is 
indicated  in  Fig.  17  by  using  water  as  the  working  fluid,  where 
the  value  of  flow  Reynolds  number  is  170.  As  can  be  seen  from 
Fig.  17,  the  present  simulation  results  agree  very  well  with  the 
experimental  and  thus  validate  the  present  model.  Their  permis¬ 
sible  difference  may  be  caused  by  the  different  position  and  size 
of  the  gas  inlets  at  the  gas  feed  headers  in  the  present  model 
and  the  experimental  by  Huang  et  al.  As  mentioned  in  Ref.  [6], 
the  validated  flow  models  can  be  established  to  simulate  gaseous 
transport  phenomena  and  electrochemical  reactions  in  the  one-cell 
stack. 

4.  Conclusions 


Dimensionless  position  along  air  flow  direction 


Fig.  16.  Averaged  temperature  difference  of  PEN  structure  for  counter  flow  in  cases 
4,  7-9  listed  in  Table  7. 


Numerical  simulation  and  modelling  the  heat  and  mass  trans¬ 
fer  and  electrochemical  reaction  in  the  one-cell  stack  of  planar 
SOFCs  are  carried  out  by  using  the  temperature-dependent  physical 
parameters  of  multi-components  mixture  gas  and  cell  stack  com¬ 
ponents.  This  study  investigates  the  degree  of  flow  uniformity  in 
different  cases  with  the  variation  of  operating  conditions  and  struc¬ 
tural  parameters  and  their  influence  to  the  cell  stack  performance. 
The  major  conclusions  are  listed  as  follows. 

(1)  Compared  with  the  co-flow  case,  the  counter  flow  case,  offers 
the  thermodynamic  and  electrochemical  advantages  in  better 
uniform  current  density  and  temperature  distributions  of  PEN 
structure  in  the  width  direction  of  the  cell  stack  due  to  the  rela¬ 
tively  uniform  air  and  fuel  flows,  as  well  as  higher  power  output, 
fuel  utilization  factor  and  fuel  efficiency. 

(2)  For  counter  flow  case,  it  is  effective  to  decrease  the  thermal 
gradient  of  PEN  structure  by  increasing  the  air  inlet  delivery 
rate  and  decreasing  the  fuel  inlet  delivery  rate.  Moreover,  better 
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performance  (higher  power  output,  fuel  utilization  factor  and 
fuel  efficiency)  is  achieved  with  lower  fuel  inlet  delivery  rate. 

(3)  For  counter  flow  case,  it  is  beneficial  to  obtain  the  better  ther¬ 
modynamic  and  electrochemical  performance  by  increasing  the 
air  inlet  temperature,  when  the  cell  stack  operates  on  the  con¬ 
dition  that  the  ascending  of  the  average  temperature  gradient 
of  PEN  structure  is  less  than  its  descending.  The  temperature 
gradient  of  PEN  structure  drops  with  the  decrease  of  fuel  inlet 
temperature. 

(4)  For  counter  flow  case,  it  is  effective  to  improve  the  thermal  gra¬ 
dient  of  PEN  structure  by  decreasing  the  anode  porosity  due 
to  the  better  uniform  air  flows.  The  power  output,  fuel  utiliza¬ 
tion  factor  and  fuel  efficiency  increase  when  anode  pore  size 
and  thickness  decrease.  But  the  temperature  distribution  of  PEN 
structure  is  almost  constant  when  anode  pore  size  decreases. 
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